classdef Economy
    % Economy Class
    properties
        state
        bank
        creditor
        government
        stats
    end
    
    methods
        function obj = Economy(objState,objBank,objCreditor,objGovernment,varargin)
            p = inputParser;
            p.addRequired('state',@(x)isa(x,'State'));
            p.addRequired('bank',@(x)isa(x,'Bank'));
            p.addRequired('creditor',@(x)isa(x,'Creditor'));
            p.addRequired('government',@(x)isa(x,'Government'));
            p.addParameter('stats',struct(),@(x)isa(x,'struct'));
            parse(p,objState,objBank,objCreditor,objGovernment,varargin{:});    

            obj.state = p.Results.state;
            obj.bank = p.Results.bank;
            obj.creditor = p.Results.creditor;
            obj.government = p.Results.government;
            obj.stats = p.Results.stats;
        end
        
        function outEcon = getEquityValue(inEcon,flagGenEq)
            
            % Inputs
            B   = inEcon.government.qBond;      % Quantity of government bonds issued at t=0
            z   = inEcon.bank.ProjPayoff;       % Payoff of project at t=2 if successful ("z")
            ell = inEcon.bank.Liquidation;      % Liquidation value at t=1 ("ell")
            LR  = inEcon.bank.LiqRatio;         % Liquidity ratio ("lambda")
            RR  = inEcon.bank.ResRatio;         % HQLA-to-deposit ratio ("rho")
            d   = inEcon.bank.qDebt;            % Quantity of debt issued by bank ("delta")
            b   = inEcon.bank.qBond;            % Quantity of government bonds held by bank ("beta")
            pE  = inEcon.bank.ProbEfficient;    % Cutoff level for the success probability of project
                                                % below which liquidation of the project is efficient("pE")

            % Partial derivatives of b and d w.r.t. lambda and rho
            den     = (LR-RR)^2/ell;
            Dd_LR   = -1/den;
            Db_LR   =-RR/den;
            Dd_RR   = 1/den;
            Db_RR   = LR/den;

            % Franchise value ("V^E", eq. 2) (= expected payoff of the 
            % risky project under the efficient continuation policy)
            [pnod,pwgt] = inEcon.state.quadState('BoundLeft',pE);
            HE = inEcon.state.cdfProb(pE);
            VE = ell*HE+z*(pwgt'*pnod);

            % Run threshold ("p*", eq. 21)
            [pR,DpR_d,DpR_b] = inEcon.bank.getProjProbAtRun();
            
            % Liquidity multiplier ("mu", eq. 31)
            LiqMult = -DpR_b/DpR_d;

            % Cost of inefficient liquidations
            if pR<1
                [pnod,pwgt] = inEcon.state.quadState('BoundLeft',pR);
                HR = inEcon.state.cdfProb(pR);
                hR = inEcon.state.pdfProb(pR);
                LiqCost = VE-ell*HR-z*(pwgt'*pnod);
            else
                HR = 1;
                hR = 0.01;
                LiqCost = VE-ell;
            end
            DLiqCost_d = (pR*z-ell)*hR*DpR_d;
            DLiqCost_b = (pR*z-ell)*hR*DpR_b;
            
            % Creditor's utility from holdings of liquid claims (eq. 13-14)
            if flagGenEq
                inEcon.creditor.LiquHoldings = B-b+d*(1-(1-LR)*HR);
            end
            LiqHold = inEcon.creditor.LiquHoldings;
            [Psi,DPsi] = inEcon.creditor.getLiquBenef();
            
            % Government bond premum ("Sb", eq. 26)
            if flagGenEq
                Sb = -(pR*z-ell)*hR*DpR_b+DPsi*(HR-(d-(ell+b))*hR*DpR_b);
            else
                Sb = inEcon.stats.Sb;
            end

            % Bank debt money premia ("Sd", eq. 27)
            Sdd     = DPsi*(d-(d-(ell+b))*HR);
            DSdd_d  = DPsi*(1-HR-(d-(ell+b))*hR*DpR_d);
            DSdd_b  = DPsi*(HR-(d-(ell+b))*hR*DpR_b);
            
            % Bank's initial equity value ("V", eq. 28)
            NetLiqBenef     = Sdd-Sb*b; 
            DNetLiqBenef_d  = DSdd_d; 
            DNetLiqBenef_b  = DSdd_b-Sb;

            V       = VE-LiqCost+NetLiqBenef;
            DV_d    = -DLiqCost_d+DNetLiqBenef_d;
            DV_b    = -DLiqCost_b+DNetLiqBenef_b;
            DV_iso  = DSdd_d*LiqMult+DSdd_b-Sb;

            DV_LR   = DV_d*Dd_LR+DV_b*Db_LR;
            DV_RR   = DV_d*Dd_RR+DV_b*Db_RR;
            
            % Welfare ("W", eq. 38)
            W = VE-LiqCost+Psi;
            
            % Output
            outEcon = inEcon;
            outEcon.stats.LiqRatio      = LR;
            outEcon.stats.ResRatio      = RR;
            outEcon.stats.qDebt         = d;
            outEcon.stats.qBond         = b;
            outEcon.stats.W             = W;
            outEcon.stats.V             = V;
            outEcon.stats.DV_d          = DV_d;
            outEcon.stats.DV_b          = DV_b;
            outEcon.stats.DV_iso        = DV_iso;
            outEcon.stats.DV_LR         = DV_LR;
            outEcon.stats.DV_RR         = DV_RR;
            outEcon.stats.LiqCost       = LiqCost;
            outEcon.stats.VE            = VE;
            outEcon.stats.NetLiqBenef   = NetLiqBenef;
            outEcon.stats.LiqMult       = LiqMult;
            outEcon.stats.HE            = HE;
            outEcon.stats.HR            = HR;
            outEcon.stats.pR            = pR;
            outEcon.stats.Sd            = Sdd/d;
            outEcon.stats.Sb            = Sb;
            outEcon.stats.LiqHold       = LiqHold;
            outEcon.stats.Psi           = Psi;
            outEcon.stats.DPsi          = DPsi;
            
        end
        
        function [FOC_RR,outEcon,outoptim] = optLiqRatio(inEcon,varargin)
            % This method solves for the bank's choice of liquidity ratio
            % given a fixed hqla-to-deposit ratio
            
            % Check if bank bond demand is bounded
            if inEcon.stats.Sb<=inEcon.stats.DPsi
                error('bank bond demand is unbounded.')
            end
            
            % Inputs
            p = inputParser;
            p.addParameter('ResRatio',double.empty,@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('TolX',10e-12,@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('TolFun',10e-12,@(x)isnumeric(x)&&isscalar(x));
            p.parse(varargin{:});
            
            ResRatio= p.Results.ResRatio;
            TolX    = p.Results.TolX;
            TolFun  = p.Results.TolFun;
            
            % Initial economy
            if ~isempty(ResRatio)
                inEcon.bank.ResRatio = ResRatio;
            end
            RR = inEcon.bank.ResRatio;
            
            % Liquidity ratio associated with 100% run probability
            LR_min = inEcon.bank.getIsoLiqRatio(1);
            
            % Optimization
            options = optimset('TolX',TolX,'TolFun',TolFun,'Display','None');
            
            fun = @(x)getFOC_LRRR(inEcon,x,RR,'LR',false);
            dx = 10e-6;
            x0 = [LR_min+dx,1-dx];
            
            funL = fun(x0(1));
            funR = fun(x0(2));
            if funR>0
                LR = x0(2);
            elseif funL<0
                LR = x0(1); 
            else
                [LR,~,~,outoptim] = fzero(fun,x0,options);
            end
            
            % Update economy with bank's choice of liquidity ratio
            [~,outEcon] = fun(LR);
            FOC_RR = outEcon.stats.DV_RR;
            
        end

            function [FOC_LR,outEcon,outoptim] = optResRatio(inEcon,varargin)
            % This method solves for the bank's choice of hqla-to-deposit ratio
            % given a fixed liquidity ratio
            
            % Check if bank bond demand is bounded
            if inEcon.stats.Sb<=inEcon.stats.DPsi
                error('bank bond demand is unbounded.')
            end

            % Inputs
            p = inputParser;
            p.addParameter('LiqRatio',double.empty,@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('x0',-3,@(x)isnumeric(x));
            p.addParameter('TolX',10e-10,@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('TolFun',10e-10,@(x)isnumeric(x)&&isscalar(x));
            p.parse(varargin{:});
            
            LiqRatio= p.Results.LiqRatio;
            x0      = p.Results.x0;
            TolX    = p.Results.TolX;
            TolFun  = p.Results.TolFun;
            
            % Initial economy
            if ~isempty(LiqRatio)
                inEcon.bank.LiqRatio = LiqRatio;
            end
            LR = inEcon.bank.LiqRatio;
            
            % Is the optimal liquidity buffer zero?
            FOC_RR0 = getFOC_LRRR(inEcon,LR,0,'RR',false);
            
            % If not, how high is the optimal liquidity buffer?
            options = optimoptions('fsolve','TolX',TolX,'TolFun',TolFun,'Display','None');
            map = @(y)LR/(1+exp(-y));
            fun = @(x)getFOC_LRRR(inEcon,LR,map(x),'RR',false);
    
            if FOC_RR0>TolFun
                [xs,~,~,outoptim] = fsolve(fun,x0,options);
            else
                xs = -inf; outoptim = struct();
            end

            % Update economy with bank's choice of hqla-to-deposit ratio
            [~,outEcon] = fun(xs);
            FOC_LR = outEcon.stats.DV_LR;
            
        end

        function [FOC_LRRR,outEcon,outoptim] = optLiqResRatio(inEcon,varargin)
            % This method solves for the bank's joint choice of 
            % liquidity ratio and hqla-to-deposit ratio
            
            % Check if bank bond demand is bounded
            if inEcon.stats.Sb<=inEcon.stats.DPsi
                error('bank bond demand is unbounded.')
            end
            
            % Inputs
            p = inputParser;
            p.addParameter('x0',-3,@(x)isnumeric(x));
            p.addParameter('TolX',10e-12,@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('TolFun',10e-12,@(x)isnumeric(x)&&isscalar(x));
            p.parse(varargin{:});
            
            x0      = p.Results.x0;
            TolX    = p.Results.TolX;
            TolFun  = p.Results.TolFun;
            
            % Is the optimal liquidity buffer zero?
            inEcon.bank.ResRatio = 0;
            [FOC_RR,outEcon,outoptim_LR] = optLiqRatio(inEcon,'TolX',TolX,'TolFun',TolFun);
            
            % If not, how high is the optimal liquidity buffer?
            options = optimoptions('fsolve','TolX',TolX,'TolFun',TolFun,'Display','None');
            map = @(y)1/(1+exp(-y));
            fun = @(x)optLiqRatio(inEcon,'ResRatio',map(x),'TolX',TolX,'TolFun',TolFun);

            if FOC_RR>0
                [xs,~,~,outoptim_RR] = fsolve(fun,x0,options);
                inEcon.bank.ResRatio = map(xs);
                [~,outEcon,outoptim_LR] = optLiqRatio(inEcon,'TolX',TolX,'TolFun',TolFun);
            else
                outoptim_RR = struct();
            end

            outoptim=cell(1,2); outoptim{1}=outoptim_LR; outoptim{2}=outoptim_RR;
            FOC_LRRR = [outEcon.stats.DV_LR,outEcon.stats.DV_RR];
            
        end
        
        function [pR,outEcon,outoptim] = optIsoRisk(inEcon,varargin)
            % This method solves for the bank's choice of liquidity ratio
            % and hqla-to-deposit ratio for a given liquidity risk appetite
            
            % Check if bank bond demand is bounded
            if inEcon.stats.Sb<=inEcon.stats.DPsi
                error('bank bond demand is unbounded.')
            end

            % Inputs
            p = inputParser;
            p.addParameter('LiqRatio',double.empty,@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('x0',[-6, 6],@(x)isnumeric(x));
            p.addParameter('TolX',10e-10,@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('TolFun',10e-10,@(x)isnumeric(x)&&isscalar(x));
            p.parse(varargin{:});
            
            LiqRatio= p.Results.LiqRatio;
            x0      = p.Results.x0;
            TolX    = p.Results.TolX;
            TolFun  = p.Results.TolFun;
            
            % Initial economy
            if ~isempty(LiqRatio)
                inEcon.bank.LiqRatio = LiqRatio;
            end
            LR = inEcon.bank.LiqRatio;
            
            % Is the optimal liquidity buffer zero?
            FOC_RR0 = getFOC_LRRR(inEcon,LR,0,'iso',false);
            
            % If not, how high is the optimal liquidity buffer?
            options = optimset('TolX',TolX,'TolFun',TolFun,'Display','None');
            map = @(y)LR/(1+exp(-y));
            fun = @(x)getFOC_LRRR(inEcon,LR,map(x),'iso',false);

            if FOC_RR0>TolFun
                [xs,~,~,outoptim] = fzero(fun,x0,options);
            else
                xs = -inf; outoptim = struct();
            end

            % Update economy with bank's choice of hqla-to-deposit ratio
            [~,outEcon] = fun(xs);
            pR = outEcon.stats.pR;
            
        end

        function [pR,outEcon,outoptim] = optLiqRisk(inEcon,varargin)
            
            % Check if bank bond demand is bounded
            if inEcon.stats.Sb<=inEcon.stats.DPsi
                error('bank bond demand is unbounded.')
            end
            
            % Inputs
            p = inputParser;
            p.addParameter('x0',[-6, 6],@(x)isnumeric(x));
            p.addParameter('TolX',10e-10,@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('TolFun',10e-10,@(x)isnumeric(x)&&isscalar(x));
            p.parse(varargin{:});

            x0      = p.Results.x0;
            TolX    = p.Results.TolX;
            TolFun  = p.Results.TolFun;
            
            % Initial economy
            pR = inEcon.stats.pR; 
            pE = inEcon.bank.ProbEfficient;

            % Optimization
            options = optimset('TolX',TolX,'TolFun',TolFun,'Display','None');
            if pR>pE
                map = @(y)1/(1+exp(-y));
                fun = @(x)optIsoRisk(inEcon,'LiqRatio',map(x),'TolX',TolX,'TolFun',TolFun);
                fun_root = @(x)(fun(x)-pR);
                [xs,~,~,outoptim_LR] = fzero(fun_root,x0,options);
            else
                outoptim_LR = struct();
                xs = inf;
            end

            % Update economy with bank's choice of liquidity ratio
            [pR,outEcon,outoptim_RR] = fun(xs);

            outoptim=cell(1,2); outoptim{1}=outoptim_LR; outoptim{2}=outoptim_RR;
            
        end
        
        function [FOC_d,outEcon,outoptim] = solveEquilibrium(inEcon,varargin)
            
            % Inputs
            p = inputParser;
            p.addParameter('x0',[-6,6],@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('TolX',10e-12,@(x)isnumeric(x)&&isscalar(x));
            p.addParameter('TolFun',10e-12,@(x)isnumeric(x)&&isscalar(x));
            p.parse(varargin{:});
            
            x0      = p.Results.x0;
            TolX    = p.Results.TolX;
            TolFun  = p.Results.TolFun;
            
            % Banks hold all the bonds in equilibrium
            ell = inEcon.bank.Liquidation;
            B   = inEcon.government.qBond;
            b   = B;
            
            % Solve for optimal quantity of debt issued by banks
            options = optimset('TolX',TolX,'TolFun',TolFun,'Display','None');
            map = @(y)ell+b+exp(y);
            fun = @(x)inEcon.getFOC_db(map(x),b,'d',true);

            [xs,~,~,outoptim] = fzero(fun,x0,options);

            % Update economy with bank's choice of debt
            [FOC_d,outEcon] = fun(xs);
            
        end
    end
    
    methods (Access = private)
        function [FOC,outEcon] = getFOC_LRRR(inEcon,LR,RR,strVar,flagGenEq)
            % This method returns the F.O.C. of the bank's choice of LR or/other RR

            inEcon.bank.LiqRatio = LR;
            inEcon.bank.ResRatio = RR;

            outEcon = getEquityValue(inEcon,flagGenEq);
            
            if strcmp(strVar,'iso')
                FOC = outEcon.stats.DV_iso;
            elseif strcmp(strVar,'LR')
                FOC = outEcon.stats.DV_LR;
            elseif strcmp(strVar,'RR')
                FOC = outEcon.stats.DV_RR;
            elseif strcmp(strVar,'LRRR')
                FOC = [outEcon.stats.DV_LR,outEcon.stats.DV_RR];
            else
                error('strVar must take the value iso, LR, RR or LRRR.')
            end
        end

        function [FOC,outEcon] = getFOC_db(inEcon,d,b,strVar,flagGenEq)
            % This method returns the F.O.C. of the bank's choice of d or/other b

            inEcon.bank.LiqRatio = (inEcon.bank.Liquidation+b)/d;
            inEcon.bank.ResRatio = b/d;

            outEcon = getEquityValue(inEcon,flagGenEq);
            
            if strcmp(strVar,'d')
                FOC = outEcon.stats.DV_d;
            elseif strcmp(strVar,'b')
                FOC = outEcon.stats.DV_b;
            elseif strcmp(strVar,'db')
                FOC = [outEcon.stats.DV_d,outEcon.stats.DV_b];
            else
                error('strVar must take the value d, b or db.')
            end
        end
    end
    
end